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ABSTRACT 

We present 3D Newtonian simulations of the coalescence of two neutron stars, using 
a Smoothed Particle Hydrodynamics (SPH) code. We begin the simulations with the two 
stars in a hard, circular binary, and have them spiral together as angular momentum is 
lost through gravitational radiation at the rate predicted by modeling the system as two 
point masses. We model the neutron stars as hard polytropes (7 = 2.4) of equal mass, and 
investigate the effect of the initial spin of the two stars on the coalescence. The process of 
coalescence, from initial contact to the formation of an axially symmetric object, takes only 
a few orbital periods. Some of the material from the two neutron stars is shed, forming a 
thick disk around the central, coalesced object. The mass of this disk is dependent on the 
initial neutron star spins; higher spin rates resulting in greater mass loss, and thus more 
massive disks. For spin rates that are most likely to be applicable to real systems, the 
central coalesced object has a mass of 2AMq, which is tantalizingly close to the maximum 
mass allowed by any neutron star equation of state for an object that is supported in part 
by rotation. Using a realistic nuclear equation of state, we estimate the temperature of 
the material after the coalescence. We find that the central object is at a temperature of 
~ lOMeV, whilst the disk is heated by shocks to a temperature of 2-4MeV. 

1. INTRODUCTION 

It has been suggested that 7-ray bursts arise at cosmological distances in the merger 
of binaries consisting of either two neutron stars, or a black hole and a neutron star 
(Eichler et. al. 1989, Piran, 1990, Paczyhski, 1991, Piran et. al. 1992, Narayan, Paczyhski 
& Piran 1992). Two different processes may provide the electromagnetic energy for the 
bursts: neutrino-antineutrino annihilation into electron-positron pairs during the merger, 
and magnetic flares generated by Parker instabilities in a postmerger differentially-rotating 
disk. Merging binary neutron stars are also a source of gravitational radiation, whose 
waveform and amplitude can be predicted with reasonable confidence (see e.g. Lincoln 
& Will, 1990; Cutler, et. al. , 1993) and are prime candidates for detection by LIGO 
(Schutz, 1986, Thorne, 1987; Abramovici et. al. 1992). It has been estimated that ~ 3 
such mergers will occur in a year within 200 Mpc (Narayan, Piran & Shemi 1991, Phinney 
1991). These events should be observable by the Laser Interferometer Gravitational Wave 
Observatory (Abramovici et. al. 1992). A merger of two neutron stars may lead to the 
ejection of neutron-rich material that may be a good site for the production of r-process 
elements. It is possible that such mergers account for all heavy r-process material in the 
Galaxy (Eichler et. al. 1989). 

The production of binaries containing two neutron stars has been discussed in the 
literature (see for examples, Narayan, Piran & Shemi 1991; van den Heuvel 1993). For 
completeness, the scenario is outlined here. Begin with two massive stars in a binary. The 
primary - the more massive of the two stars - will evolve first, forming a neutron star 
in a supernova. If the binary is tight enough, and the secondary sufficiently massive, the 
binary may not be disrupted by this violent event. Mass transfer between the secondary 
and neutron star will ensue as the former evolves beyond the main sequence, the star 
either filling its Roche lobe, or losing mass via a wind. At this stage in its evolution, the 
system will manifest itself as a massive X-ray binary. If the secondary were to evolve to 
a supernova without significant mass loss and if the explosion is symmetric, the binary 



would be assured of disruption, as the mass lost from the exploding star would be much 
greater than that remaining in the two neutron stars. However, we may speculate upon the 
following panacea: because the bloated secondary is more massive than the neutron star, 
mass transfer from the former will be unstable, if the mass transfer rate is large enough, 
bringing the two stars closer together until they form a common envelope system. Such a 
system represents an efficient way to remove the nefarious excess matter in the secondary's 
envelope. By finely tuning the separation of the neutron star and core of the secondary 
at the beginning of the common envelope phase, we are able to speculate that one may 
remove a large fraction of the gaseous envelope, leaving the core of the secondary and the 
neutron star in a tight binary. When the core finally explodes, the relative mass loss from 
the system may thus have been reduced sufficiently to avoid disruption. Alternatively, 
the supernova explosion might be assymetric leaving the neutron stars in an eccentric 
orbit. It is unclear how likely the binary is to remain intact. Estimates based on the 
current known neutron stars binaries suggest that this happens in about one percent of 
the cases (Narayan, Piran & Shemi, 1991). The fragility of such systems is evidenced by 
the eccentricities of the four known neutron star binaries, suggesting these systems were 
close to disruption when they were formed. It should be noted that systems formed in this 
manner will have circularized before the two neutron stars come into contact (C. Cutler, 
private communication) . 

Two neutron stars in a binary will spiral together as angular momentum is lost via 
the emission of gravitational radiation. The timescale for merger of two neutron stars, 
with Ml — M2 = IAMq, is less than 10^° years if the initial separation satisfies d ^ 5Rq 
(Iben and Tutukov 1984). Recent work (Kochanek 1992; Bildsten & Cutler, 1992; Lai, 
Rasio & Shapiro, 1993) has shown that the viscosity of the neutron-star fluid is unlikely to 
be sufficient to spin-up the stars to achieve tidal-locking as the two stars spiral together. 
One therefore has to ask how the merger process will change as a function of neutron star 
spin. 

Hydrodynamic simulations of neutron-star mergers are a formidable task. Beyond 
the 3D hydrodynamic calculations it should include general relativitistic effects, employ a 
realistic equation of state, contain neutrino transport and neutrino cooling, and possible 
nuclear reactions. So far there have been several attempts to address various aspects of 
the problem. Most of these attempts focused on the gravitational radiation emmission 
(Oohara & Nakamura, 1989, 1990; Nakamura & Oohara 1989, 1991; Shibata, Nakamura 
& Oohara, 1992). These were Newtonian calculations using a finite difference code, a 
polytropic equation of state and including a gravitational radiation backreaction formula. 
Rasio & Shapiro (1992) focused on the hydrodynamic evolution of the rotating core and 
began with the two neutron stars corotating at contact. Here we present a study that 
focuses on the thermodynamics and nuclear physics of the coalescence. As was mentioned 
above, the low viscosity of the neutron-star material leads us to consider the merger of non- 
corotating neutron stars. In order to achieve a realistic estimate of the radial velocity at 
the moment that the neutron star come in to contact we begin with the two neutron stars 
being far apart and follow the spiral-in phase. We follow the evolution until the neutron 
stars merge and we examine both the matter distribution and the thermal properties of the 
system after the two neutron stars merge. Our goal is to estimate the thermal conditions 
in the coalesced core and the disk that forms around it after the two neutron stars merge. 



and to study the possible subsequent nuclear processes and neutrino radiation processes. 
As a first stage we do not attempt to include relativistic effects, or even the complete 
equation of state in these calculations. At this stage, we perform Newtonian calculations 
with a polytropic equation of state, with a gravitational radiation backreaction formula. 
We employ, however, a realistic equation of state to estimate the thermodynamic conditions 
during the computations. We compare the conditions prevailing in our simulations to those 
required to produce a gamma-ray burst by neutrino-antineutrino annihilation (Eichler et. 
al. 1989) or by accretion of the disk onto the compact core (Narayan, Paczyhski & Piran 
1992). 

The numerical methods applied to this problem and the initial conditions used are 
described in section 2. The results of our simulations are presented in section 3 and in 
section 4 we discuss the various physical processes that could take place in this system. 
The implications of these results to the question of formation of gamma-ray bursts and to 
r-process nucleosynthesis are summarized in section 5. 



2. NUMERICAL METHODS AND INITIAL CONDITIONS 



We performed a series of simulations of the merger of two neutron stars of equal mass, 
using a 3D Smoothed Particle Hydrodynamics (SPH) code, based on one used previously to 
study collisions between stars (see for example, Benz & Hills 1992). Being a Langrangian 
particle code, SPH is well suited to this problem. We have no need for a computational box 
(gravitational forces being computed using a tree), hence we do not waste computational 
resources simulating the evolution of the voids between the neutron stars as they spiral 
together. 

Since we are mostly interested in the thermodynamic conditions of the material we are 
using a Newtonian code, rather than a general relativistic one. The only general relativistic 
effect which we include is the addition of the gravitational radiation backreaction force, 
which we discuss below. We modelled the cold neutron stars as hard polytropes. We 
employed a realistic nuclear equation of state developed by Lattimer & Swesty (1991), to 
compute pressure as a function of density for cold material at nuclear densities with a 
compressibility K = 180 MeV, and obtained 7 = 2.4 for a polytropic approximation of the 
equation of state. We checked aposteriori that the Lattimer & Swesty equation of state is 
well approximated by a polytrope for the density range that appears in our calculations 
(except for a small amount of mass in the very low density region). In view of the large 
uncertainty in the equation of state of nuclear matter we believe that this is a reasonable 
first step. We later make use of the Lattimer & Swesty equation of state to estimate the 
temperature from the internal energy and density at the final stages of the computation. 
To construct the initial neutron-star model, we solve the Lane-Emden equation for 7 = 2.4, 
and we construct a polytrope with the required density profile using the method described 
in Davies, Benz & Hills (1992). 

We began each simulation with the two stars 5i?ns apart and had them spiral in at 
the rate predicted by treating the system as two point masses and considering the energy 
lost through the emission of gravitational radiation. The change in energy, E, and angular 



momentum, J, is given by (e.g. Shapiro & Teukolsky 1983) 

(1) 
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where a is the binary separation, M the total mass of the system, and fj, = MiM2/{Mi + 
M2) is the reduced mass, where Mi = M2 = (1/2)M is the mass of the two neutron stars. 
Using the above equations, we derive the foUowing formulae for the acceleration due to the 
emission of gravitational radiation. 
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As we apply the same acceleration to each SPH particle within each neutron star, the 
circulation of the fluid will be conserved, as the accelerating field has a zero curl. This 
in-spiral force was turned off once the two stars came into contact. At this stage it is 
clear that this approximation breaks down. While this force is a far cry from a sofisticated 
gravitational radiation emission backreaction force (Blanchet, Daniour & Schafer 1990) it 
provides us with a simple formula that can be easily added to our Newtonian calculations 
which will provide a realistic radial velocity component at the moment when the two stars 
collide. 

As the two stars spiral together, a lag angle will develop between the stars' long-axes 
and the line joining their center s-of-mass. The size of this angle is partially a function of 
the viscosity of the neutron star fluid. As the fluid inside the neutron stars is extremely 
inviscid, we expect the flgures of the two neutron stars to remain essentially alligned, until 
the very last stages of the in-spiral. However, one might be concerned that the much larger 
effective viscosity present in our fluid would lead to an artiflcially-large lag angle, with the 
two stars then receiving torques. Indeed, in early, lower-resolution runs, the stars were 
spun up close to the orbital spin rates. To circumvent this problem, we ran a second set of 
mergers, alligning the two neutron stars every time-step, so that their long-axes remained 
parallel. The shift in total angular momentum due to this procedure was measured and 
found to be tolerably small (~ 1 in 10^). In fact, the two higher resolution runs presented 
here (runs C and D in Table I) that used this reallignment mechanism show very similar 
results to the equivalent runs performed without the reallignment (runs A and B). In other 
words, the problem which seemed to require the use of reallignment largely disappeared 
when the resolution was increased in the flnal runs we performed. Both sets of runs are 
given here for completeness. 

Each simulation was run for 3.5ms (~ 50 dynamical times) after the stars came into 
contact. This was sufficient time for all mergers to produce a central coalesced object 



surrounded by a disk. As we show later, our simplifying assumptions do not break down 
at this stage. However, the evolution of this system on longer time scales is determined by 
processes, which we did not include, but address in detail in section 4. 

3. SPH POLYTROPIC RESULTS 

We performed four high-resolution mergers (runs A - D) and two lower resolution 
simulations (runs E and F). In runs A - D, we used 4271 particles per neutron star, in 
runs E and F, we used 1256 particles per neutron star. In these six runs we keep all initial 
parameters the same except the initial spins of the two neutron stars. The spin angular 
velocities of the neutron stars when they came into contact is given in Table 1. For run 
A the lag angle ~ 4° by the time the two stars come into contact. For run B, the value 
is somewhat larger at ~ 10°. Taking Mns = 1.4M0 and R^s = 10km, Ogpin = 1 (the 
synchronous spin with the two neutron stars 2i?ns apart) corresponds to a spin period of 
0.46ms. Hence in runs A and C, the two neutron stars have spin periods ~ 2.5ms. We 
consider this to be a reasonable upper limit to the neutron star spins for any merging 
system, as the stars are unlikely to be spun up to their orbital periods just before contact, 
owing to the low viscosity of the neutron star fluid. In runs B and D, we consider the 
case where the two neutron stars have negligible spin. We expect essentially all systems 
to lie in between these two cases. In runs E and F, we further investigate the effect of the 
neutron star spins, by considering cases where the stars are spinning in the opposite sense 
to the orbital spin. 

3.1 Phenomenology of Mergers 

In all cases we observe the appearance of a single rotating coalesced object surrounded 
by a disk. The size of the disk varies greatly from one run to another. In some cases 
additional long tidal tails form as well. Density contour plots for runs C and D are given 
in Figures 1 and 2. Logarithmic contours are plotted at intervals of 0.25 dex, beginning 
at a density of 0.001 in code units (1 code unit = 2.786 x 10^^ g/cm"^). The sequences 
shown begin just after the onset of coalescence, times being given in code units (1 time 
unit = 73/Us) with t = when the two stars were 5i?ns apart. In both cases, a single 
coalesced object forms on the timescale of one orbit {i.e. ~ 1ms) after the objects collide. 
The system quickly becomes axisymmetric after the two stars come into contact, hence 
only a small amount of energy will be emitted as gravitational radiation at this stage, as 
has been calculated in earlier simulations by Shibata, Nakamura & Oohara (1992) and 
by Rasio & Shapiro (1992). Hence our approximation of switching off the gravitational 
radiation when the two stars come into contact has been vindicated. 

In run C (where the neutron stars are rapidly rotating), far larger tidal tails develop 
behind the neutron stars. These tails are clearly visible by t = 306, and are less visible 
at the same time for run D. By a time, t = 314, some of the material removed from the 
neutron stars in the tidal tails has begun to form a disk around the central object. As 
we will see shortly, this material is shock-heated. By the end of the simulation, a large 
fraction of the mass removed from the two stars remains in the disk: only some of the 
material thrown off initially remains, unshocked, in the two tidal tails. We thus observe a 
difference between runs C and D: in the latter case, the lower neutron-star spin results in 
a delay in the formation of tidal tails, the material has thus been ejected less far when the 
two objects coalesce, and any material that has been removed from the two neutron stars 
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is located in the disk by the end of our simulation. 

Figure 3 illustrates the distribution of SPH particles, in the plane of the original 
orbits, at the end of Run C. This plot illustrates the three components of the final system 
produced in our simulation: the central coalesced object (out to a radius ~ 2Rns), a disk 
(out to a radius ~ QRns)^ and two elongated tidal tails, which have expanded to ~ 20i?ns 
by the end of the simulation. One has to address different questions regarding each of these 
regions. We wish to know whether the energetics of the core or the disk are sufficient to 
fuel cosmological gamma-ray bursts. For this we need to examine on one hand the thermal 
energies and the temperatures of the central coalesed object and the disk as well as the 
mass of the disk, and on the other hand the opacity for neutrino emission in various regions. 
We are concerned with the stability and subsequent evolution of the central object and for 
this we are interested in the rotational and gravitational energies of the core. Finally we 
must also estimate how much material is ejected from the system. This question has some 
implication to the emission of gamma-ray bursts by these events and to the possiblity that 
r-process nucleosynthesis is taking place in these sites. 

In Figure 4, we plot density contours of the final configuration for run C in a plane 
perpendicular to that of the neutron star trajectories. Logarithmic contours are plotted, 
at intervals of 0.25 dex, starting at a density of 0.0001 in code units (recalling 1 code 
unit = 2.786 x 10^^ g/cm^). The maximal density decreases by about 10% as the two 
stars are tidally distorted and then subsequently increases by about 5% relative to the 
initial density as a central coalesced object forms. The final central density is ~ 5 times 
the nuclear density. This is only slightly higher than the central density of the original 
neutron stars. At these densities the polytropic index derived from the Lattimer & Swesty 
(1991) equation of state is practically unchanged. Figure 4 shows the form of the central, 
coalesced object and the disk of material surrounding it. It is clear from this plot that the 
disk is thick, almost toroidal; the material having expanded on heating through shocks. 
This disk surrounds a central object that is somewhat flattened due to its rapid rotation. 
It is apparent from this flgure that an almost empty centrifugal funnel forms around the 
rotating axis and there is almost no material above the polar caps. The funnel is also 
apparent in Figure 5, where we produce a contour plot of the integrated column density, 
i.e. J^pdz. This funnel might be important for the formation of gamma-ray bursts as 
it provides a region in which a baryon free radiation-electron-positron plasma could form. 
Such plasma is an essential step in any cosmological gamma-ray burst model. 

3.2 Angular Velocity and Mass Distribution 

For each run, we computed the ratio of system kinetic energy to the gravitational 
potential energy, i.e. T/\W\, as a, function of time. All runs show somewhat similar 
behavior, with the initially more spun-up neutron stars having slightly greater kinetic 
energy. In all runs, T/|VF| increased from ~ 0.1 to ~ 0.17 between t = 290 and t = 300. 
As mass was shed from the neutron stars when they produced tails, T/|VF| decreased, 
reaching a value of ~ 0.13 at the end of the simulations. For rigidly-rotating polytropes 
of 7 ^ 2.2, a secular instability occurs for T/|M^| = 0.14. Under such an instability, a 
rotationally-symmetric ellipsoid (Maclaurin spheroid) will transform into a triaxial ellipsoid 
(Jacobi ellipsoid). However, such a transformation will occur on a relatively slow timescale 
(Press and Teukolsky 1973), i.e. many dynamical timescales. Hence, even though the 
polytrope considered in our simulations is sufficiently hard, and the system reaches a large 



enough value of T/|M^|, the swirhng mass has resolved itself into a central coalesced object 
surrounded by a thick disk before the system has time to transform into a Jacobi ellipsoid. 
A discussion concerning the subsequent stability of the central coalesced object is given in 
section 4. 

In Figure 6 we plot the spin angular velocities of the SPH particles as a function of 
cylindrical radius. For both runs C and D, again, we see evidence for the three regions 
produced in run C. The material in the central object (r ^ 1.5i?ns) rotates rigidly; an 
effect of the viscosity of our SPH fluid (which is much larger than that of the real neutron 
star fluid). In run D, the rotation curve is somewhat broadened in the same region. This is 
due to the lower spin rates of the material within the neutron stars before they coalesced. 
It is unlikely that the real neutron star fluid will be as orderly within the coalesced object, 
as it is so inviscid. For 1.5i?ns ~ t ^ 6Rns {i-e. the disk), Ogpin decreases almost as r~^/^; 
the slight flattening of the power law caused by the mass contained in the disk. For run C, 
the material located at greater radii is contained in the two tails. Note that, as discussed 
earlier, no tails remain in run D, and a mere few, stray particles, are found at such large 
distances from the center-of-mass. The material in the tails produced in run C closely 
follows a power law, O, oc r~^, where (3 ~ 1.8, i.e. the rotational velocity falls off faster 
than Keplerian. This may, initially, seem a rather enigmatic result, for which we furnish 
the following explanation. The essential point to realize is that the particles contained in 
the tails are travelling on eccentric paths. Imagine flrst the simple case where all these 
particles were ejected with the same angular velocities and at the same distance from the 
center-of-mass, but at different times. We would then have an ensemble of particles located 
at different positions on the same ellipse - or rather a set of identical ellipses with different 
orientations of the semi-major axes in the orbital plane. Kepler's Second law tells us that 
Or^ = a constant for a given ellipse, hence if our simple model were correct, we would 
expect to see O oc r~^. However, imagine now that the particles ejected flrst (and now 
located at the largest radii) were ejected with more angular momentum. Their "constant" 
would be therefore larger than that for particles ejected later, and thus now closer to the 
central object. This will have the effect of softening the power law, as is observed. 

In Figure 7 we plot the enclosed mass (as a fraction of total mass) as a function of 
cylindrical radius for all runs. Again, we see the three components of the systems produced 
in Runs A and C: the central object contains ~ 85% of the system mass, the disk some 
13%, with the remaining 2% of the material being found in the elongated tidal tails at 
Tcyi ^ GRns- Runs B and D show somewhat similar proflles to those of A and C, except 
they leave no material in tidal tails, and the disk is slightly less massive, containing ~ 10% 
of the system mass. For the two cases where the neutron stars were spinning initially in 
the direction opposite to the orbital motion, Mdisk ~ O.OSMsystem- The dependence of this 
function on the initial neutron star spins is clear: higher spin rates (in the same sense as 
orbital spin) result in greater mass-loss, and more massive disks. 

3.3 Temperatures and Thermal Energy 

In Figure 8 we plot the energy per unit mass of each SPH particle against the 
density as the particle site for run C (both are in code units: units of energy/mass being 
1.858 X 10^° ergs g~^). We computed isotherms using Lattimer and Swesty's nuclear 
equation of state, assuming nuclear matter in beta equilibrium, and were thus able to 
estimate the temperatures at the SPH-particle locations. As stated earlier, this is an 




approximate method, which is particularly sensitive for the high density regions in which 
the isotherms are very near each other, that is a small change in u produces a large change 
in T. The isotherms derived from this method are shown in Figure 9. A lot of the material 
around the Z-axis in the final configuration has come from outer regions of the neutron 
stars, has thus been compressed considerably, and thus heated, increasing its T. Typical 
temperatures in the central, coalesced object are around T ~ 10 — 15MeV. It should be 
noted that even if we are overestimating the temperature of the central coalesced object 
by a factor of two here, the conclusions drawn regarding the timescale for neutrino escape 
in section 4.1 are not radically altered. When material is tidally shredded from the two 
neutron stars, it initially decompresses adiabatically, then subsequently shock heats as the 
material from the two tidal tails come into contact, forming a disk. This shocked material 
now has densities in the range 3 x lO^^g cm~"^ ^ p ^ 3 x 10^"^g cm~^. Again the same 
polytropic index for the EOS applies as for higher densities, with the exception of the 
very outermost low density layers. In the simulation, the disk is heated by shocks to a 
temperature of 2-4 MeV, and contains ~ 0.25 x 10^^ ergs of kinetic energy, and 2 x 10^^ 
ergs in thermal energy. 

For reasonable uncertainties in the equation of state, Lattimer and Yahil (1989) 
obtained a relation between gravitational mass and binding energy of non-rotating neutron 
stars, which reads 



Et^n = (1.5 ± 0.15) f ^ ) X lO^Vg. (3) 



If we apply it to two 1.4 Mq neutron stars and a 2.4 Mq central object (assuming it is 
stable), we find a difference of roughly 2.76 xlO^^ erg in binding energy which would be 
expected to go into internal energy (heat), as long as the system did not loose energy by 
radiation yet. However, the system is not in its ground state, but in an excited state due to 
its rotation. The binding energy is decreased by the rotation energy. In the Newtonian limit 
we have Erot = 0.5/0^, with / denoting the moment of inertia. In an order of magnitude 
estimate we make use of the moment of inertia for a homogeneous object {OAMR^ for a 
sphere and O.bMR^ for a disk). The central coalesced object is slightly elongated and we 
take a medium factor 0.45, a radius of 17km (cutting at p = 5 x 10^"^g cm~"^, see Figure 
4), and fl — 8.6 x 10"^s~^. This leads to i?rot=2.32xl0^"^erg or a reduced binding energy 
increase between two single neutron stars and the rotating coalesced object of 4.4xl0^^erg 
or 9.2 xlO^^ erg g~^ = 9.55 MeV per nucleon, which should show up as additional thermal 
energy in the formed compact object and is only due to stronger gravitational binding. 
This agrees well with the temperature (internal energy) of the central object of 10-15 MeV 
(see Figure 9). 

This excercise also illustrates that the object is a relatively fast rotator (although not 
at break-up and rotationally supported) , which leads to the fact that its central density is 
only increased by roughly 5% over that of the initial neutron stars (to about 1.7 < lO^^g 
cm~^). In that regime the same polytropic index applies to the Lattimer and Swesty (1991) 
EOS as at lower densities and the object is still stable. 
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4. ADDITIONAL PHYSICAL PROCESSES 

The previous sections discussed the hydrodynamical behavior during the first few 
miUiseconds of the merger of two neutron stars, resulting in one central object with densi- 
ties up to 10^^ g cm"^ and a disk, containing about O.2-O.4M0 of matter with densities of 
3x10^^ — 1 X lO^^g cm~"^. This part of the calculation included only (adiabatic) hydrody- 
namics and a partial loss of gravitational radiation. There are two questions which have 
to be asked at this point: (i) do energy loss or energy production mechanisms exist which 
work on shorter time scales, thus making the present result invalid, and (ii) which phys- 
ical processes will shape the future behavior of that system, its stability, energy balance 
between nuclear energy generation and neutrino losses, and its composition. We address 
these questions below. 

4.1 Neutrino Escape 

Let us first assume that all constituents are in thermal equlibrium and that 2 MeV 
neutrinos populate the disk and 10 MeV neutrinos the central object. The neutrino nucleon 
scattering cross section (see e.g. Tubbs and Schramm 1975, Shapiro and Teukolsky 1983) 
is given by 



1 / E.^ 



2 



cr.,n = ^c7o{^^) ^ao^^(MeV) ctq = 1.76 x lO'^^m^ (4) 

V A/ tpO 



The mean free path between scattering events Xi, = l/{na^^n) depends on the scat- 
tering cross section a^^n and the nucleon number density n — pNa- This can be expressed 
in the following form 



nn ,^5 /lOMeVV lO^^gcm" ,, 

X^ = 9.9 X 10^ cm. (5) 

\ E„ J p 



Thus, for a typical density of 10^^ g cm~"^ and 2 MeV neutrinos in the disk we find 
a mean free path of 246 km and for an average density of the central object of about 
8 X 10^^ g cm~^ and 10 MeV neutrinos we find 12.4 m. In case the disk settled already to 
a nuclear statistical equilibrium for that temperature and density, it will consist partially 
of alpha particles (see Lamb et al. 1978), which leads however to the same result (we have 
to multiply Xj^ by ~ N'^/A = 4/ A, making use of the coherent scattering cross section for 
nuclei and approximating the Weinberg angle by sin^ 6*^^=0. 25). 

When scattering events follow a random walk in three dimensions, the timescale for 
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travelling an absolute distance d is 

t{d) = ^. (6) 

It has been pointed out by Burrows and Lattimer (1986) and Burrows (1988) that neutrinos 
do not follow a scattering random walk, but that their behavior is more complex and 
absorptions can occur as well. Nevertheless, equation (6) gives a good order of magnitude 
estimate, sufficient for our discussion. When employing a typical height of the disk of 15 
km and 10 km for the central object, this leads to typical diffusion time scales of 9.1 x 10~^s 
in the disk and 8.1 x 10~^s in the central object. The first number indicates that neutrino 
escape would already be substantial in the disk at a time of 1ms, while this is not the case 
for the central object that cools on a longer time scale. The difference with respect to a 
supernova core collapse, where typical neutrino diffusion time scales of about Is occur, is 
due to the somewhat higher temperatures obtained in the proto neutron star after stellar 
collapse and bounce. With thermal energies being higher by about a factor of 2.5-3.0, 
leading to a factor of 6-9 in mean free path and diffusion time scale, where the energies 
enter quadratically, we have tdiff = 0.5 — 0.7s. 

4.2 Neutrino Emission and Cooling 

Now we have to consider which processes produce these neutrinos and whether their 
production time scales are as fast as the loss time scales. The very high densities and 
high temperatures of the central object produce neutrinos of all families essentially instan- 
taneously in comparison to the diffusion time scales. Neutrinos exist in thermal equilib- 
rium and all neutrino families carry similar energies. The total internal energy gain of 
^ 4.4 X lO^^erg will be released with a diffusion time scale of roughly 0.08s (as derived 
above). This differs from a supernova core collapse, where the proto neutron star releases 
about 2.5 X lO^'^erg on time scales of Is. 

Accretion onto the proto neutron star in a supernova core collapse with a rate M, 
before the formation of the delayed shock via neutrino heating, will increase the neutrino 
luminosity by (see e.g. Burrows 1988) 




1 + ^1 - 2GM/c^Ry ,^. 

M \ f M \ /I0km\ _i 

erg s 



O.IMqs-i; \MqJ \ R J 



A similar accretion phase, fed by matter in the disk, can occur onto the central object 
in systems of merged neutron star binaries. 

While in the central object energy loss due to neutrinos is limited by the diffusion 
time scale controlling the possible escape, we noticed that the neutrino diffusion timescale 
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in the disk is very short (of the order of 10~^s). All neutrino poduction processes occur 
on longer time scales, and thus a free escape for neutrinos from these neutrino cooling 
processes is guaranteed. Positron capture on neutrons and electron capture on protons 
which emit anti-neutrinos and neutrinos dominate over all other neutrinos losses (pair, 
photo, plasma, and bremsstrahlung neutrinos - Schinder et al. 1987, Itoh et al. 1989, 
1990). The rates for the disk conditions are of the order of 10^*^ — lO^^erg g~^s~^ (Fuller et 
al 1980,1982). The next largest energy loss rate due to pair and plasma neutrinos amounts 
only to about 10^^ — lO^^erg g~^s~^ (Itoh et al. 1989, 1990). The time scale for electron and 
positron captures (and accompanying neutrino production) is of the order 10~^ — 10~^s, 
which underlines our previous conclusion of free escape (with diffusion time scales of the 
order of 10~^s), i.e. the full energy loss is encountered and no trapping occurs. 

4.3 Thermonuclear Energy Generation in the Disk 

In addition to the cooling processes, we have to consider nuclear energy generation 
in the disk. Simple nuclear network calculations show that any initial protons combine 
with neutrons to form alpha particles on a time scale of less than 10~^^s for the density 
and temperature conditions in the disk. This releases an energy of 7.07 MeV per nucleon 
or 6.8 X lO^^erg g~^ or 1.36 x 10^^ erg per O.IM© of disk matter burned to helium. Disk 
material in beta equilibrium would consist of less than 1% protons. Thus, a disk in 
weak equilibrium would burn 2% of its mass (equal amounts of neutrons and protons) 
instantaneously, releasing about lO^'^erg (on a 10~^^s time scale). This would heat and 
somewhat expand the disk because of a higher energy generation than neutrino cooling 
rate. 

The subsequent evolution is complex and depends on many details, i.e. the resulting 
density and temperature of the disk and the free proton fraction Yp. High densities (large 
electron Fermi energies) and/or temperatures favor electron captures on protons. For the 
disk densities electron capture time scales of about 10~^ — 10~'^s with a neutrino cooling 
rate close to lO^'^ — lO^^erg g~^s~^ are obtained in the temperature range of 1-lOMeV. The 
positron capture rates on neutrons are generally lower, because positrons have a negative 
Fermi energy. Similar values as for electron captures on protons can only be attained at 
high temperatures and lower densities. These rates would usually balance in order to keep 
beta equilibrium. But each proton produced will burn to helium on much shorter time 
scales, thus dominating over electron captures. The depletion in the free proton abundance 
Yp is replenished by further positron captures on neutrons, and gradually the free neutrons 
would burn to He via conversion to protons. 

A total energy gain would occur if the density and proton abundance (determining 
the electron Fermi energy and electron capture rate on protons) become small enough that 
neutrino losses from electron captures are negligible, and the temperatures are low enough 
that the neutrino losses in each positron capture reaction (producing one proton) become 
smaller than 14MeV, which is 1/2 of the energy released when two resulting protons burn 
with two available neutrons to He. This is the case for temperatures of <lMeV, pY^ < lO^g 
cm~^, and Yp < 10""^ (tables 2 and 3 in Fuller et al. 1982 and private communication). 
The low temperatures are required as the weak capture cross sections scale with E'^ and 
thus the average energy of the positron and the escaping neutrino are higher than the 
thermal value. For the conditions discussed above the time scale of burning neutrons via 
protons to He is of the order 100s. 
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On the other hand the low densities ensure electron Fermi energies of less than 5MeV. 
Thus, the formation of heavier elements via neutron captures and beta decays with typical 
Q-values of 13-15MeV will not be inhibited by electron captures on neutron-rich nuclei. An 
r-process similar to the r-process suggested for inhomogeneous big bang nucleosynthesis 
(Applegate, Hogan, Scherrer 1988; Thielemann et al. 1991) with light nuclei embedded 
in a neutron bath will occur. Then the long time scales required to convert neutrons 
via protons to helium can be overcome by this feature of the burning process. Once 
the first heavy seeds are available, neutron captures and beta decays (not temperature 
dependent) will drive the burning on beta decay time scales of 0.1 to 0.01s and cause a net 
energy deposition in the low density and temperature parts of the disk. The maximum 
energy of photons emitted in the nuclear reactions would result from beta decays far from 
stability with decay energies up to 15-20 MeV. The scenario is very similar to that of "cold 
decompression" (Lattimer et al. 1977, Meyer 1989), with the exception that no energy loss 
terms and an adiabatic expansion fueled by the energy generation was considered in the 
latter. The total expected energy release is of the order / x (5 — 6) x 10^^ erg, dependent 
on the distribution of nuclei produced and the fraction / of the disk involved. 

In the interior of the colliding neutron stars one finds 5-30% protons, depending 
on the equation of state (Weber and Weigel 1989). This increase is due to the fact that 
neutrons and protons are also degenerate at such high densities and the neutron and 
proton Fermi energies are of importance as well. If such matter would be ejected during 
the collision on shorter than weak equilibrium time scales, up to 60% of the disk matter 
could burn instantanesously, releasing about 3 x lO^^erg in the disk. Our calulations 
show that the ejection of matter into the disk occurs on time scales comparable to weak 
equilibrium time scales and a behavior closer to the earlier discussion is expected. 

In none of the cases would the energy release be sufficient to unbind and "evaporate" 
the whole disk, which has a total binding energy of ~ 4 x lO^^erg but possibly the outer 
lower density parts can escape. Further calculations which treat energy loss rates and 
thermonuclear reactions consistently as part of the hydro code will hopefully be able to 
answer the open "details" . 

4.4 Stability of the Central Coalesced Object 

Our hydro calculations produced a central object which was stable on the time scales 
considered in this calculation. We have to consider several points to understand that 
outcome and its generality. The main question concerns the limiting neutron star mass 
for a given equation of state. The object is rotating slower than, but close to the, Kepler 
frequency of typically Qk = 10^s~^. This allows an increase beyond the limiting mass of 
non-rotating neutron stars of 14-17% (Friedman, Ipser, Parker 1986; Friedman and Ipser 
1987; Weber and Glendenning 1992), the largest values are related to the softest equations 
of state. 

The core is rapidly rotating but Jc/GM^ ^ .6 < 1. (similar results were obtained by 
Shibata Nakamura & Oohara (1992) and by Rasio & Shapiro (1992)). Thus the core is not 
supported completely by rotation. This can be understood intuitively from the following 
argument. The angular momentum just before the stars are in contact is 

J « 2MOi?2 = 2M{GM/4R^)^/^R^ = (GM^R)^/^. (8) 
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The ratio J / {AGM'^ / c) determines whether the object is rotationally supported. We find 



Jc U R V" 



AGW^ 4 \rg{M) 



(9) 



where rg{M) is the gravitational radius of one of the stars, i.e. r^ ^ 2{M/lAMQ)kui. With 
R f» lOkm we have that R/vg ~ 5 and hence Jc/GM^ ^ .5. The ratio is sUghtly higher 
in our case because of the spin of the neutrons stars. However, it is quite impossible to 
reach Jc/GM^ > 1. The intuitive explanation for the phenomenon is that while orbiting 
each star is attracted only by its companion and the centrifugal force has to balance this 
gravitational attraction. However, once the two objects coalesce the gravitational force 
towards the center results from a total mass of 2M which is too large to be stopped purely 
by the rotation. 

The behavior of the equation of state at high densities is dominated by the com- 
pressibility K and the possible appearance of new degrees of freedom, i.e. existence of 
excited baryons (hyperons) in addition to neutrons and protons, kaon condensates, etc. If 
it is permitted to populate these states rather than neutrons and protons with their huge 
Fermi energies, a smaller pressure and softer EOS results, which causes a smaller limiting 
neutron-star mass. Our polytropic index was adjusted to the Lattimer and Swesty EOS 
with a compressibility K = 180 MeV. Only neutrons, protons, and leptons are considered 
at high densities. This corresponds to a large limiting neutron star mass of 2.2M0. 

Most recent relativistic equations of state result in limiting masses for non-rotating 
neutron stars below 2.2M0, down to 1.4-1. GM© (Glendenning et al. 1992). With a 14-17% 
increase for objects close to their Kepler frequencies, we would start out with stable central 
objects for the upper end of this mass range. We would expect neutrino release from the 
hot object and mass accretion until accretion, cooling or loss of angular momentum lead to 
surpassing the limiting mass and cause the collapse to a black hole. We expect that this will 
take place on a time scale of seconds. A limiting non-rotating neutron star mass smaller 
than 2M0 would lead to an unstable coalesced object from the time of its formation. Thus, 
a black hole would form after about 1ms and a neutrino burst of only a few times lO^'^erg 
would result during the collapse (Gourgoulhon and Haensel 1993). This will be enhanced 
by continuing accretion from the disk. 

This clearly shows that the outcome of a neutron star merger, its neutrino signal, 
electromagnetic radiation, and probably also the disk behavior depend decisively on the 
limiting neutron star mass. If identified with observed events (i.e. a gamma ray burst), we 
could finally expect observational contraints for this theoretically quite uncertain quantity. 

5. DISCUSSION 

For all the neutron-star spins considered here, the merger results in a rapidly rotating, 
though not completely rotationally supported, central coalesced core surrounded by a disk 
of ejected material. A few percent of the total mass is ejected into long extended tails 
in some of the cases. During the merger event material from the two neutron stars is 
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compressed adiabatically, increasing the temperature of the gas to ~10 MeV and it becomes 
a strong source of neutrinos. The neutrino coohng time is around 10 — 100ms. The hot 
rapidly rotating core is on the edge of stabiUty and its subsequent evolution depends on 
the true equation of state. The disk material gets heated via shocks and its temperature is 
^ 5MeV. Since its density is much lower its neutrino cooling time is significantly shorter. 
Thermonuclear reactions can, however, heat the disk and its thermal evolution depends 
on a delicate balance between heating and cooling. 

Our calculations included several simplifying approximations and it is worth while to 
examine their validity before turning to the interpretation of the results: (i) We employed a 
polytropic equation of state instead of a realistic one. However we found that by the end of 
the calculations the central density increased only by 5% and the polytropic approximation 
to the Lattimer & Swesty (1991) equation of state is valid at this density range, (ii) We 
did not include neutrino transport and cooling. The neutrino cooling time of the core 
(10 — 100ms) is sufficiently long that there are no significant neutrino losses up to the 
moment that we stop the calculations. The time scale for neutrino cooling of the disk 
is, however, much shorter and our approximation breaks down there, (iii) We switch 
off the gravitational radiation backreaction force when the two neutron stars touch each 
other. The core becomes axisymmetric within a fraction of a millisecond after the two 
stars have come in contact, and only a small amount of energy is released via gravitational 
radiation during this phase (Shibata et. al., 1993, Rasio & Shapiro 1993). (iv) Our worst 
approximation is that apart from gravitational radiation backreaction we did not include 
other general relativistic effects. There is no physical justification for that, only a technical 
one, the required three dimensional hydrodynamics calculations are beyond the current 
scope of numerical relativity. General relativistic effects are the major uncertainty in the 
calculations that we present here. We have no way of estimating them quantitatively. 
Qualitatively we expect that they might enhance the dynamical instability of the central 
rotating core. 

The mass of the central object ranges between 2.4 — 2.QMq. This is above most 
upper limits for cold neutron star masses even with the stiffest equations of state, (e.g. the 
limiting mass is 2.2M0 for non-rotating neutron star with the Lattimer- Swesty equation 
of state that we follow). It is clear that if the real equation of state is soft, the central 
object must be unstable and it will collapse quickly to a black hole. The situation is less 
clear with a hard equation of state. The maximum mass allowed increases by up to 17% 
when the neutron star is rotating rapidly. Furthermore, our hydrodynamics calculations 
show that the Newtonian system with a hard equation of state, seems stable and does 
not show any tendency to collapse. At the end of the calculation the central polytropic 
index is large which indicates that the system is far from the instability that sets in as 
the polytropic index decreases to 4/3. Here it is possible that the combination of rotation 
and thermal pressure may prevent the direct collapse. The collapse might be induced later 
when additional mass from the disk accretes onto the central object or when it cools down. 
General relativistic effects, which we ignore here, might have a critical effect and cause the 
central merged object to collapse sooner. As the central core is on the edge of stability 
it might be the case that in some cases (depending e.g. on the spin angular momentum) 
the core remains extended initially while in others it collapses directly to a black hole. In 
either case the ultimate final configuration is a black hole. 
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The total neutrino emission depends on how fast the core collapses to a black hole. 
One can expect appreciable neutrino emission if the rotating configuration is stable, or 
if the combination of rotational and thermal energy make the system stable, halts the 
collapse until the core cools or accretes additional matter. If general-relativistic effects or 
a softer equation of state will make the system unstable a black hole will form after a few 
msec without much neutrino emission form the core. 

The disk extends up to 60km or so. Its total mass is ~ OAMq, and its total kinetic 
energy ^ 2.5 x 10^^ erg. From our SPH runs, we see that the disk contains ~ 10^^ erg 
of thermal energy, with some of the material being shock-heated up to temperatures of ~ 
2-4 MeV. The disk is optically thin to neutrinos, and neutrino cooling exceeds heating for 
temperatures T >lMeV, densities pY^ > lO^g cm~^ and free proton fraction Yp > 10~^. 
Hence we expect neutrino cooling to lead to a fiattened disk in the inner high density part. 

There occurs a thermonuclear fiash of burning of the initial protons (with neutrons) in 
the disk to He. At locations in the disk where after that initial fiash T, pYf, ,and Yp are equal 
or smaller than the values listed in the previous paragraph net energy deposition rather 
than neutrino losses will dominate and a further expansion is expected. The low densities 
ensure electron Fermi energies less than 5MeV. Thus, the build-up of heavy elements from 
the seed nuclei, produced in the early thermonuclear fiash, via neutron captures and beta 
decays (with Q-values of typically 13-15MeV) is not inhibited by reverse electron captures 
on these neutron-rich nuclei. If a fraction / of the whole disk is burned to heavy elements, 
this would result in a total release of / x (5 — 6) x 10^^ erg. We expect some of the material 
of the outer disk and extended arms to become unbound from the system due to this 
nuclear energy release and to contain r-process matter produced during the decompression 
phase. 

Ratnatunga & van den Bergh (1989) find core collapse supernova rate in our galaxy 
of 2.2 ± 1.1 per century, 10~^ Mq to 10~^ Mq of r-process material must be ejected per 
supernova event (Cowan, Thielemann & Truran 1991) to produce the observed abundance 
of r-process material. Eichler et. a/.(1989) suggested that, alternatively, the r-process 
takes place in neutron star mergers. They estimated the production rate using the merger 
rate estimated by Clark van den Heuvel & Sutantyo. Recently, Narayan et. a/.(1991) and 
Phinney (1991) estimate that there are io~^-^='=-^ neutron star mergers per year per galaxy. 
Using the updated rate we find that to explain r-process nucleosynthesis as a by product 
of neutron star merger one requires the ejection of 10~'^Mq to IMq per event. The upper 
limit is clearly impossible, however, the lower range is of the same order of magnitude of 
what we observed in the long tails in some of the mergers. Some of the r-process material 
is observed in very low metalicity stars and galactic evolution models would require that 
it is produced as early as a few x 10^ — lO^years after the formation of the galaxy. This 
is comparable with the life time of the famous binary pulsar PSR-1913+16 and one might 
expect that there is no problem caused by this time scale. This suggests that some r-process 
material could be produced via this mechanism even on such a short time scale. 

We turn now to the implications to 7-ray bursts models. Cosmological 7-ray bursts 
require ~ 10^^ ergs with a rise time as short as tens of msec (for some of the short bursts). 
The hot coalesced core is one possible source of this energy. If the core does not collapse 
directly to a black hole it will emit its thermal energy as neutrinos. The neutrino fiux is 
sufficiently large that ~ 10~^ to 10""^ of it could be converted to electron-positron pairs 
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via vv -^ e+e" (Goodman et. ai, 1987) and those could produce a 7-ray burst (Eichler 
et. al., 1989). The time scale for the neutrino burst is short enough to accommodate even 
the shortest rise times observed. Additional neutrino luminosity could arise on a longer 
time scale from accretion of the disk material on the central object. 

An additional energy source that could power a 7-ray burst from a neutron star 
merger is the disk surrounding the central object. This energy source can operate regardless 
of the question of whether the central object collapse directly to a black hole or not. In 
this case the neutrino emission is not sufficient to power a burst via neutrino annihilation, 
however, Narayan et. a/. (1992) have suggested that the magnetic energy density in the 
disk can reach equipartition with the kinetic energy (corresponding to 10^^ gauss) and 
recombination of such a magnetic field could produce an electromagnetic burst. It is 
not clear what will be the time scale for the emission from this disk. It could take several 
seconds or even longer, depending on the viscosity in the disk, which is unknown at present. 
As was stressed by Narayan et. a/. (1992) the large variability in observed properties of 
7-ray bursts could be explained by those two alternative energy sources. An additional 
source of diversity, which we discover here, is the distinction between systems that collapse 
directly to a black hole and those that undergo a longer rotating core phase. 

The initial 7 rays that are produced are not observed directly. The huge optical 
depth of the resulting electron-positron plasma would result in an optically thick fireball 
(Goodman, 1986; Shemi & Piran, 1990; Paczynski, 1990; Piran & Shemi, 1993). The 
observed 7-rays emerge much later when the fireball becomes optically thin (due to its 
rapid expansion) or when it interacts with interstellar material (Meszaros & Rees, 1993). 
This late stage determines both the duration and the spectrum of the bursts. Shemi & 
Piran (1990) have shown that to produce a 7-ray burst the fireball must be practically free 
from baryonic load. If more than 1O~^M0 baryons are injected into the fireball they will 
prevent the formation of a 7-ray burst. Larger amounts of mass are ejected in a neutron 
star merger. However all the mass is ejected into the equatorial plane. In all mergers 
simulated, we note the the regions above the poles of the coalesced object are relatively 
free of material. The accuracy of our simulation is of course not enough to estimate the 
amount of matter in those funnels, but it seems that those funnels may provide an avenue 
for neutrino escape and the production of a clean radiation fireball. This suggests that the 
7-rays are beamed and appear only in certain directions. 
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Table I 
Characteristics of Runs 



Run 


' 'contact 


Npart 


Reallign 


Afdisk/Afsystem 


A 


0.2329 


8542 


no 


0.13 


B 


0.0478 


8542 


no 


0.10 


C 


0.2143 


8542 


yes 


0.13 


D 


0.0109 


8542 


yes 


0.10 


E 


-0.0439 


2512 


yes 


0.05 


F 


-0.1257 


2512 


yes 


0.05 
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FIGURE CAPTIONS 

Figure 1: Density contour plots for run C, in ttie plane of the neutron star trajecto- 
ries. Time and distances are given in code units (1 time unit = 73//s, 1 distance unit 
= Rns =10km). Logarithmic contours are plotted, at intervals of 0.25 dex, beginning at a 
density of 0.001 in code units (1 code unit = 2.786 x lO^^g/cm"^). 

Figure 2: Density contour plots for run D. Units as given in Figure 1. 

Figure 3: The distribution of SPH particles, in the plane of the original orbits, at the 
end of run C. Distance units as given in Figure 1. 

Figure 4: Density contour plot of the final configuration in run C in a plane perpendicular 
to that of the neutron star trajectories. Logarithmic contours are plotted, at intervals of 
0.25 dex, starting at a density of 0.0001 in code units (1 code unit = 2.786 x 10^^ g/cui^). 

Figure 5: Contour plot of column densities, i.e. at a given location (r, z), column density, 

Pc = /^ pdz. 

Figure 6: The spin angular velocities of the SPH particles as a function of cylindrical 
radius. For run C (figure a) and run D (figure b). Both angular velocity and distance 
are in code units (1 distance unit = Rns =10km, O = 1 is equivalent to a spin period of 
0.461ms). 

Figure 7: The enclosed mass (as a fraction of total mass) as a function of cylindrical 
radius (in units of the neutron-star radius) . 

Figure 8: The energy per unit mass of each SPH particle against the density as the 
particle site for the final configuration in run C. Both are given in code units, energy/mass 
units being 1.858 x 10^*^ ergs/gm. 

Figure 9: Isotherms produced using Lattimer and Swesty's nuclear equation of state, 
assuming beta equilibrium. Isotherms are draw for T =1.75, 2, 3, 5, 7, 10, 15, and 25MeV. 
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